## Loading required package: maptools
## Loading required package: sp
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
## Loading required package: RColorBrewer
## Loading required package: MASS
## Loading required package: rgeos
## rgeos version: 0.5-9, (SVN revision 684)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.5-0
## Polygon checking: TRUE
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
##
## Attaching package: 'raster'
## The following object is masked from 'package:MASS':
##
## select
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/user/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/user/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.4-0
##
## Attaching package: 'spatstat.geom'
## The following objects are masked from 'package:raster':
##
## area, rotate, shift
## The following object is masked from 'package:MASS':
##
## area
## Loading required package: spatstat.random
## spatstat.random 2.2-0
## Loading required package: spatstat.core
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
##
## getData
## Loading required package: rpart
## spatstat.core 2.4-4
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-2
##
## spatstat 2.3-4 (nickname: 'Watch this space')
## For an introduction to spatstat, type 'beginner'
planning_area_sf <- st_read("data/Planning_Area_Census2010.kml")
## Reading layer `Planning_Area_Census2010' from data source
## `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\Planning_Area_Census2010.kml'
## using driver `KML'
## Simple feature collection with 55 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 103.6054 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
class(planning_area_sf)
## [1] "sf" "data.frame"
Next we read in point data. The first set of point data is the Park Connector (PCN) access points themselves. Next, we have data on various public amenities like supermarkets, hawker centres, pharmacies as well as transportation infrastructure like MRT exits.
access_points_sf <- st_read("data/PCN_Access_Points.kml")
## Reading layer `PCN_Access_Points' from data source
## `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\PCN_Access_Points.kml'
## using driver `KML'
## Simple feature collection with 306 features and 2 fields
## Geometry type: POINT
## Dimension: XYZ
## Bounding box: xmin: 103.6933 ymin: 1.272978 xmax: 104.0048 ymax: 1.46272
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
supermarkets_sf <- st_read("data/supermarkets.kml")
## Reading layer `SUPERMARKETS' from data source
## `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\supermarkets.kml'
## using driver `KML'
## Simple feature collection with 742 features and 2 fields
## Geometry type: POINT
## Dimension: XYZ
## Bounding box: xmin: 103.6258 ymin: 1.24715 xmax: 104.0042 ymax: 1.462167
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
hawkercentres_sf <- st_read("data/hawkercentre.kml")
## Reading layer `HAWKERCENTRE' from data source
## `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\hawkercentre.kml'
## using driver `KML'
## Simple feature collection with 125 features and 2 fields
## Geometry type: POINT
## Dimension: XYZ
## Bounding box: xmin: 103.6974 ymin: 1.272716 xmax: 103.9882 ymax: 1.449217
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
mrt_exits_sf <- st_read("data/lta-mrt-station-exit-kml.kml")
## Reading layer `MRT_EXITS' from data source
## `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\lta-mrt-station-exit-kml.kml'
## using driver `KML'
## Simple feature collection with 474 features and 2 fields
## Geometry type: POINT
## Dimension: XYZ
## Bounding box: xmin: 103.6368 ymin: 1.264972 xmax: 103.9893 ymax: 1.449157
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
pharmacy_sf <- st_read("data/registered_pharmacy.kml")
## Reading layer `REGISTERED_PHARMACY' from data source
## `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\registered_pharmacy.kml'
## using driver `KML'
## Simple feature collection with 269 features and 2 fields
## Geometry type: POINT
## Dimension: XYZ
## Bounding box: xmin: 103.699 ymin: 1.2563 xmax: 104.0025 ymax: 1.448227
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
class(access_points_sf)
## [1] "sf" "data.frame"
class(supermarkets_sf)
## [1] "sf" "data.frame"
class(hawkercentres_sf)
## [1] "sf" "data.frame"
class(mrt_exits_sf)
## [1] "sf" "data.frame"
class(pharmacy_sf)
## [1] "sf" "data.frame"
Next, we read in raster data, which consists of only population density data.
GDALinfo("data/popmap15adj.tif")
## Warning in GDALinfo("data/popmap15adj.tif"): statistics not supported by this
## driver
## rows 368
## columns 542
## bands 1
## lower left origin.x 103.6383
## lower left origin.y 1.164905
## res.x 0.0008333
## res.y 0.0008333
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=longlat +datum=WGS84 +no_defs
## file data/popmap15adj.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float32 TRUE -3.402823e+38 128 128
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 -4294967295 4294967295 NA NA
## Metadata:
## AREA_OR_POINT=Area
pop_density <- raster("data/popmap15adj.tif")
pop_density
## class : RasterLayer
## dimensions : 368, 542, 199456 (nrow, ncol, ncell)
## resolution : 0.0008333, 0.0008333 (x, y)
## extent : 103.6383, 104.09, 1.164905, 1.471559 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : popmap15adj.tif
## names : popmap15adj
First, we plot the PCN access points. On first glance, it is apparent that these points cannot be randomly distributed, with noticeable clusters in the Choa Chu Kang/Bukit Batok area as well as the Woodlands/Sembawang area. The points also seem to occur as a “string” of points, which obviously stands to reason since these are access points to Park Connectors which are essentially lines.
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red")+ tm_layout(title="PCN Access Points in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid
Now, we will display the geographical distribution of the various amenities (supermarkets, hawker centres, etc.) on top of the locations of the PCN acces points. ### Supermarkets We plot the locations of supermarkets overlaid with the access points on a map of Singapore. We see that most access points are located very close to supermarkets.
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red") + tm_shape(supermarkets_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "skyblue") + tm_layout(title="PCN Access Points in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid
For hawker centres, the correlation is not as obvious.
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red") + tm_shape(hawkercentres_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "green") + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid
For MRT exits, we notice that the Downtown Core area has a high concentration of exits. This could be because many of the stations in that area have many exits, sometimes as many as 7 or 8. This makes it seem more numerous. We can try getting a dataset of just MRT stations themselves to help us narrow down our investigation. Unlike the MRT exits, we see that there are very few PCN access points in the Downtown Core area. This could be because of the lower priority assigned to cycling paths over mass rapid transit in a space-scarce area like the Downtown Core.
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red") + tm_shape(mrt_exits_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "orange") + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid
For pharmacies, there doesn’t seem to be a clear pattern.
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red") + tm_shape(pharmacy_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "purple") + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid
Interestingly, areas with the highest population density (Serangoon, Hougang and the areas surrounding Paya Lebar) seem to be relatively underserved by acces points as compared to Western regions like Bukit Batok which has a much lower population density.
Standing out as well is the central part of Singapore, which has virtually zero population density. This region encompasses the MacRitchie Reservoir area and naturally is not home to any human population density. It might be useful to exclude this region in later analysis to prevent confounding.
tmap_mode("view")
## tmap mode set to interactive viewing
# tm_shape(pop_density) + tm_raster() + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "green") + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
tm_shape(pop_density) + tm_raster() + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "green") + tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid
training_sf <- planning_area_sf[planning_area_sf$Name %in% c("YISHUN","JURONG EAST", "JURONG WEST", "TAMPINES", "PASIR RIS", "BEDOK"),]
training_sf
## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 103.6749 ymin: 1.296503 xmax: 103.9848 ymax: 1.445616
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
## Name
## 1 PASIR RIS
## 8 JURONG EAST
## 14 JURONG WEST
## 23 YISHUN
## 35 TAMPINES
## 52 BEDOK
## Description
## 1 <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>PASIR RIS</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>PASIR RIS</td> </tr> <tr> <td>Planning Area Code</td> <td>PR</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>EAST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>ER</td> </tr> <tr> <td>INC_CRC</td> <td>77F6CE717F532D1A</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:19 PM</td> </tr> <tr> <td>X_ADDR</td> <td>40795.475</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>40066.2774</td> </tr> <tr> <td>SHAPE_Length</td> <td>22444.999809</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>15810712.053754</td> </tr> </table> </td> </tr> </table> </body> </html>
## 8 <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>JURONG EAST</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>JURONG EAST</td> </tr> <tr> <td>Planning Area Code</td> <td>JE</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>WEST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>WR</td> </tr> <tr> <td>INC_CRC</td> <td>E0C3B6AB25F8AFFF</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:19 PM</td> </tr> <tr> <td>X_ADDR</td> <td>17035.1876</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>33633.557</td> </tr> <tr> <td>SHAPE_Length</td> <td>27342.347259</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>17857762.682791</td> </tr> </table> </td> </tr> </table> </body> </html>
## 14 <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>JURONG WEST</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>JURONG WEST</td> </tr> <tr> <td>Planning Area Code</td> <td>JW</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>WEST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>WR</td> </tr> <tr> <td>INC_CRC</td> <td>3302EAE21AD1E998</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:20 PM</td> </tr> <tr> <td>X_ADDR</td> <td>13705.4978</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>35973.1233</td> </tr> <tr> <td>SHAPE_Length</td> <td>19356.680643</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>14688294.312111</td> </tr> </table> </td> </tr> </table> </body> </html>
## 23 <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>YISHUN</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>YISHUN</td> </tr> <tr> <td>Planning Area Code</td> <td>YS</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>NORTH REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>NR</td> </tr> <tr> <td>INC_CRC</td> <td>8D4108597B11C0DC</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:19 PM</td> </tr> <tr> <td>X_ADDR</td> <td>28472.3369</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>44137.8596</td> </tr> <tr> <td>SHAPE_Length</td> <td>20778.98769</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>21578985.748052</td> </tr> </table> </td> </tr> </table> </body> </html>
## 35 <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>TAMPINES</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>TAMPINES</td> </tr> <tr> <td>Planning Area Code</td> <td>TM</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>EAST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>ER</td> </tr> <tr> <td>INC_CRC</td> <td>C5A6D638AA4E65FE</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:20 PM</td> </tr> <tr> <td>X_ADDR</td> <td>41454.9098</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>36208.5982</td> </tr> <tr> <td>SHAPE_Length</td> <td>22402.760996</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>20913646.028123</td> </tr> </table> </td> </tr> </table> </body> </html>
## 52 <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>BEDOK</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>BEDOK</td> </tr> <tr> <td>Planning Area Code</td> <td>BD</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>EAST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>ER</td> </tr> <tr> <td>INC_CRC</td> <td>DBC1F8A651648374</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:21 PM</td> </tr> <tr> <td>X_ADDR</td> <td>38582.8088</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>34032.1699</td> </tr> <tr> <td>SHAPE_Length</td> <td>21853.967373</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>21731365.846805</td> </tr> </table> </td> </tr> </table> </body> </html>
## geometry
## 1 MULTIPOLYGON Z (((103.9532 ...
## 8 MULTIPOLYGON Z (((103.7119 ...
## 14 MULTIPOLYGON Z (((103.7281 ...
## 23 MULTIPOLYGON Z (((103.8468 ...
## 35 MULTIPOLYGON Z (((103.9643 ...
## 52 MULTIPOLYGON Z (((103.9636 ...
sg_area <- 728.6 #km2
supermarkets_gd <- 742 / sg_area
hawkercentres_gd <- 125 / sg_area
mrt_exits_gd <- 474 / sg_area
pharmacy_gd <- 269 / sg_area
sg <- readOGR("./data", "Region_Census2010")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data", layer: "Region_Census2010"
## with 5 features
## It has 8 fields
sg
## class : SpatialPolygonsDataFrame
## features : 5
## extent : 2637.869, 56396.44, 15748.72, 50256.33 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs
## variables : 8
## names : OBJECTID, REGION_N, INC_CRC, FMEL_UPD_D, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area
## min values : 1, CENTRAL REGION, 51D899B637F87BAF, 2014/05/06, 13007.354, 31931.7133, 60174.6538079, 112912795.479
## max values : 5, WEST REGION, C05F166F6716C12E, 2014/05/06, 42244.19, 44177.9547, 238296.73891, 251603415.909
sg_window <- as.owin(sg)
# convert crs of access_points_sf to crs of sg
access_points_q <- st_transform(access_points_sf, crs(sg))
# convert sf to ppp
access_points_ppp <- as.ppp(st_coordinates(access_points_q), W=sg_window)
## Warning: data contain duplicated points
Try plotting quadrats of various sizes.
Q6 <- quadratcount(access_points_ppp, nx=6, ny=6)
plot(Q6)
Q10 <- quadratcount(access_points_ppp, nx=10, ny=10)
Q15 <- quadratcount(access_points_ppp, nx=15, ny=15)
Q20 <- quadratcount(access_points_ppp, nx=20, ny=20)
Q40 <- quadratcount(access_points_ppp, nx=40, ny=40)
plot(Q40)
From the quadrat density graph, we again see that most access points are
concentrated in the North and West.
# compute density within each quadrat
Q_density <- intensity(Q40)
# plot the density raster
plot(intensity(Q10, image=TRUE), main=NULL, las=1)
# plot access points
plot(access_points_ppp, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)
# plot the density raster
plot(intensity(Q20, image=TRUE), main=NULL, las=1)
# plot access points
plot(access_points_ppp, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)
# plot the density raster
plot(intensity(Q40, image=TRUE), main=NULL, las=1)
# plot access points
plot(access_points_ppp, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)
pop_density_q <- projectRaster(pop_density, crs = crs(sg))
brks <- c(0, 50, 100, 150, 200, 250, 300)
zcut <- cut(pop_density_q, breaks = brks)
pop_tess <- tess(image = zcut)
Q2 <- quadratcount(access_points_ppp, tess = pop_tess)
## Warning in quadratcount.ppp(access_points_ppp, tess = pop_tess): Tessellation
## does not contain all the points of X
plot(Q2, main="", las=1)
# compute density within each quadrat
Q2_density <- intensity(Q2)
Q2_density
## tile
## 1 2 3 4 5 6
## 3.029003e-07 1.125329e-06 1.760026e-06 5.016187e-07 4.759687e-07 5.626379e-07
cl <- interp.colours(c("blue", "green", "yellow", "orange", "red", "darkred"), 6)
# plot the density raster
plot(intensity(Q2, image=TRUE), main=NULL, las=1, col=cl)
# plot access points
plot(access_points_ppp, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)
# Monte Carlo Q20 (Shu Ling)
ann_observed <- mean(nndist(access_points_ppp, k=1))
ann_observed
## [1] 300.9152
n <- 1000 # number of simulations
ann_pop <- vector(length = n) # empty object to be used to store simulated ANN values
pop_density_im <- as.im(pop_density_q) # convert population density to im object
for (i in 1:n) {
random_point <- rpoint(n=access_points_ppp$n, f=pop_density_im)
ann_pop[i] <- mean(nndist(random_point, k=1))
}
Window(random_point) <- sg_window
plot(random_point, pch=16, main=NULL, cols=rgb(0,0,0,0.5))
hist(ann_pop, main=NULL, las=1, breaks=40, col="bisque", xlim=range(ann_observed, ann_pop))
abline(v=ann_observed, col="blue")
## MRT exits
# convert crs of access_points_sf to crs of sg
mrt_exits_q <- st_transform(mrt_exits_sf, crs(sg))
# convert sf to ppp
mrt_exits_ppp <- as.ppp(st_coordinates(mrt_exits_q), W=sg_window)
MRT_Q20 <- quadratcount(mrt_exits_ppp, nx=20, ny=20)
plot(MRT_Q20)
MRT_density20 <- intensity(MRT_Q20, image=TRUE)
plot(MRT_density20)
class(MRT_density20)
## [1] "im"
n <- 1000 # number of simulations
ann_mrt <- vector(length = n) # empty object to be used to store simulated ANN values
for (i in 1:n) {
random_point <- rpoint(n=access_points_ppp$n, f=MRT_density20)
ann_mrt[i] <- mean(nndist(random_point, k=1))
}
Window(random_point) <- sg_window
plot(random_point, pch=16, main=NULL, cols=rgb(0,0,0,0.5))
hist(ann_mrt, main=NULL, las=1, breaks=40, col="bisque", xlim=range(ann_observed, ann_mrt))
abline(v=ann_observed, col="blue")